EPJ manuscript No. 

(will be inserted by the editor) 



OO 
O 

o 

(N 

c 

00 
(N 



o 
u 

d 



> 
in 

(N 
(N 

o 

oo 
o 



Capillary filling for multicomponent fluid using the 
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Abstract. We present a systematic study of capillary filling for a binary fluid by using mesoscopic a lattice 
Boltzmann model describing a diffusive interface moving at a given contact angle with respect to the 
walls. We compare the numerical results at changing the ratio the typical size of the capillary, H, and 
the wettability of walls. Numerical results yield quantitative agreement with the theoretical Washburn 
law, provided that the channel height is sufficiently larger than the interface width and variations of the 
dynamic contact angle with the capillary number are taken into account. 

PACS. 83.50.Rp , - 68.03.Cd 



1 Introduction 

The physics of capillary filling is an old problem, originat- 
ing with the pioneering works of Washburn T] and Lucas 
j2]- It remains, however, an important subject of research 
for its relevance to microphysics and nanophysics 4,5, 
[3]. Capillary filling is a typical "contact line" problem, 
where the subtle non-hydrodynamic effects taking place 
at the contact point between liquid-gas and solid phase 
allow the interface to move, pulled by capillary forces and 
contrasted by viscous forces. As already remarked, Wash- 
burn in 1921 [1\ described theoretically the dynamics of 
capillary rise. Considering also inertial effects, except the 
"vena contracta" , and two fluids with the same density 
(p a = 1, pb = p a ) and the same viscosity (jj, a = p b = A*), 
the equation of motion of the moving front is [3] : 



d 2 z(t) 12 p, dz(t) _ 2-y cos6 



dt 2 



H 2 p dt 



HpL 



(1) 



used is a suitable adaptation of the Shan-Chen pseudo- 
potential LBE [7] with hydrophobic/hydrophilic bound- 
aries conditions, as developed in [5]. 



2 LBE for capillary filling 

The relevant geometry is depicted in fig. fl}. The bot- 
tom and top surfaces are coated only in the right half of 
the channel with a boundary condition imposing a given 
static contact angle |S]; in the left half we impose peri- 
odic boundary conditions at top and bottom surfaces in 
order to mimic an "infinite reservoir" . Periodic boundary 
conditions are also imposed at the two lateral sides such 
as to ensure total mass conservation inside the system. At 
the solid surface, bounce back boundary conditions for the 
particle distributions were imposed. The conditions which 
allow the wetting of the surfaces will be discussed in the 
following. 



where H is the capillary height, L its length, 7 the surface 
tension and 9 the contact angle. This model is obtained 
under the assumption that (i) the instantaneous bulk pro- 
file is given by the Poiseuille flow, (ii) the microscopic slip 
mechanism which allows the motion of the interface is not 
relevant to bulk quantities (such as the overall position 
of the interface inside the channel), (iii) inlet and out- 
let phenomena can be neglected (limit of infinitely long 
channels). In the following, we will show to which extent 
this phenomenon can described by a mesoscopic Lattice- 
Boltzmann equation for multicomponent. The model here 



2.1 LBE algorithm for multi-component flows 

Let us review the multicomponent LB model proposed 
by Shan and Chen [7]. This model allows for distribution 
functions of an arbitrary number of components with dif- 
ferent molecular mass: 



At 
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f?(x+ Ci At,t+At)-tf(x,t) = /*(*,*) - f* (eq) (x,t) 

Tk L 

where ff(x,t) is the kinetic probability density function 
associated with a mesoscopic velocity Cj for the fcth fluid, 
Tfc is a mean collision time of the fcth component (with At 
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a time lapse), and f^ eq \x,t) the corresponding equilib- 
rium function. For a two-dimensional 9-speed LB model 
(D2Q9) fi {eq \x,t) takes the following form [9]: 

rk(eq) 2 e „ e „ 

/ = a k n k - -n k u k H ■ u k H 

Meq) (1 - a k )n k 1 eq 

h - — 5 — + z nkCi ' u k 

+ ^n k (ci ■ u e k q f - ^n k u e k q ■ u e k q for i=l. . .4 (3) 
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n k c t ■ \i k H 



-n k (a ■ u q f - ^n k u e k q ■ n eq for i=5. 



(4) 



In the above equations c^'s are discrete velocities, defined 
as follows 



0,i = 0, 

/ (i-lW • (»-l)?r\ • 
c, = < [cos^-,sin^-j ,1 



= 1-4 



V2 ( C0S [^ + f ], 8in[U=fZ + | ]) 



i = 5 — 8 



( 5 ) 

in the above, is a free parameter related to the sound 
speed of a region of pure fcth component according to 
( c s) 2 = K 1 ~ a k)\ n k = Y,ifi is the number density 
of the fcth component. The mass density is defined as 
pk = m kn k , and the fluid velocity of the fcth fluid 
is defined through p k u k — m k J2i c ifi> where m k is the 
molecular mass of the fcth component. The equilibrium 
velocity u^ 9 is determined by the relation 



PkU-t 9 = Pfcu' + T k F k 



(6) 



where u' is the common velocity of the two components. 
To conserve momentum at each collision in the absence of 
interaction (i.e. in the case of = 0) u' has to satisfy 
the relation 



u 




(7) 



The interaction force between particles is the sum of a 
bulk and a wall components. The bulk force is given by 

s 

F lfe (x) = ]T J2 G kk W k (x')(x' - x) (8) 

x' fe=l 

where G kk is symmetric and tf'fe is a function of n k . In our 
model, the interaction-matrix is given by 



g kk , \x' - x\ = 1, _ 
Gkk = { g kk /±,\x'-x\ = A 
0, otherwise. 



(9) 



where g kk is the strength of the interparticle potential 
between components fc and fc. In this study, the effective 
number density •Pfe(nfc) is taken simply as tf'fc(rife) = n k . 



ny=H 




Fig. 1. Geometrical set-up of the numerical LBE. The two di- 
mensional geometry, with length 2L and width H , is divided in 
two parts. The left part has top and bottom periodic bound- 
ary conditions such as to support a perfectly flat gas-liquid 
interface, mimicking a "infinite reservoir". In the right half, of 
length L, there is the true capillary: the top and bottom bound- 
ary conditions are those of a solid wall, with a given contact 
angle 9. Periodic boundary conditions are also imposed at the 
west and east sides. 



Other choices would lead to a different equation of state 
(see below). 

At the fluid/solid interface, the wall is regarded as a 
phase with constant number density. The interaction force 
between the fluid and wall is described as 



F 2fc (x) 



-nk(x)^2 9kv 



,(x')(x' 



(10) 



where n w is the number density of the wall and g kw is the 
interaction strength between component fc and the wall. 
By adjusting g kw and n w , different wettabilities can be 
obtained. This approach allows the definition of a static 
contact angle 0, by introducing a suitable value for the 
wall density n w [5], which can span the range 8 £ [0° : 
180°]. In particular, we have chosen g% w = 0,<?2tu = —512 
while n w is varied in order to adjust the wettability. In the 
sequel, we choose 512 = 0.2 which indicates that species 
2 is attracted by the wall (hydrophilic) , while species 1 is 
neutral. Let us note that high values of associated 
with hydrophilicity. 

In a region of pure fcth component, the pressure is given 
byPfc = (c^) 2 m k n k , where (c^) 2 = |(1 — ccfe). To simulate 
a multiple component fluid with different densities, we let 
(c k s ) 2 m k = Cq, where Cg = 1/3. Then, the pressure of the 
whole fluid is given by p = c 2 , J2 k n k + § Efc.fc 9k~k^k^k, 
which represents a non-ideal gas law. The viscosity is given 
by v — §Q7Jfc PkT k — \), where (3 k is the mass density 
concentration of the fcth component. 

The Chapman-Enskog expansion [5] shows that the 
fluid mixture follows the Navier-Stokes equations for a sin- 
gle fluid: 

dtp + V • (pa) = 0, (11) 
p[d t u + (u ■ V)m] = -VP + F + V • (vp(Vu + uV) 
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Fig. 2. Front displacement vs time for different channel height 
H — 15, 30, 50, fOO with their corresponding analytical solu- 
tions. The discrepancy from washburn's law is stronger for the 
smallest channel. The channel length is always L = 450 except 
for H=f00, for which L = 500. 



Fig. 3. Front dynamics for different n w (1.0, 0.75, 0.5) that is 
for different degree of wettability. The configuration considered 
is with H = 30. The = 0.5 is fitted by an analytical 

solution with 9 = 0.78 and the case n w = 0.75 with 6 — 0.52. 



where p = J2k Pk ^ s the total density of the fluid mix- 
ture and the whole fluid velocity u is defined by pu = 

2.2 Numerical Results 

All simulations were performed using the Shan-Chen model 
described above, setting vi = v g = 0.167, pi = p g = 1, 
gi2 — 0.2, a = 4/9, that is c 2 = |, and the interfacial 
tension is 7 — 0.07. The channel length is chosen to be 
L = 450. By taking 9 constant in time, a simple analyti- 
cal solution of equation ([1]) can be obtained: 

= VcapHcos f ^ [ exp (_ t /^) + t / td _ !] + ZQj (12 ) 

where zq is the starting point of the interface at the be- 
off 2 

ginning of the simulation, td — ^77 ^ s a typical transient 
time and V cap = — is the capillary speed. This solution 
has been used to compare with simulations. 

The front displacement as a function of time is shown 
in Fig. [2] for different values of the channel height H = 
15,30,70,100, at n w = 1, for which static contact angle 
was found to be 9 ss 5. As expected, the velocity of the 
front grows with channel height. The analytical curves are 
given by the solution of Eq. (fl"2|) . where the contact angle 
is the dynamic one computed from numerical data. The 
contact angles computed for the four heights 15, 30, 50, 100 
are respectively 0°, 11°, 25°, 45°. The dynamic contact an- 
gle has been obtained directly as the slope of the contours 
of near-wall density field, and independently through the 
Laplace's law, AP = 2 ^ c ° se , The latter has been chosen 
for the comparison with analytical fitting curves, because 
the direct computation from density contours turns out 
to be less precise. Nevertheless, the values calculated in 
the two ways are approximately consistent. For instance, 
the contact angle computed for the case H = 30 from the 



direct measurements of the pressure is 9 « 12° against 
the value 9 w 11° computed via density contours. Some 
comments on the front dynamics are in order.. The case of 
smallest channel height does not follow the analytical so- 
lution, showing the finite size of the interface (w/H 1/3) 
significantly affects the results. On the other hand, for 
a larger channel, good agreement between numerical and 
theoretical results not only holds asymptotically, but it 
also extents to the initial transient. This is particularly 
true for the largest height H = 100, where the transient 
time-scale td = is sufficiently long to make the expo- 
nential term in the solution (|12|) important over a macro- 
scopic time span. The results show that the dynamic con- 
tact angles experience a strong dependence on the channel 
height. In particular, in small channels, dynamic contact 
angles remain near their static values. On the other hand, 
for large ones the discrepancy is evident. This is due to 
the increasing value of the capillary number (Ca ~ 0.03 
for H = 100), since it is known that there is a correc- 
tion of the dynamic contact angles due to finite capillary 
numbers. This correction takes the form the general form 
cos(9d) — cos(9 s ) = g(Ca). Our results are best fitted by 
g(Ca) — 18 Ca 12 , which is in line with previous forms 
used in different LBE methods [TT1IT2"] 

Hereafter the configuration with H — 30 and n w — 
1.0 will be used as a reference for all simulations. In fig- 
ure [H the front dynamics is shown for the 
1.0, 0.75, 0.5. As expected, it is found that more hy- 
drophobic cases correspond to smaller velocities . The ana- 
lytical solutions which fit the numerical data are obtained 
respectively with 9 = 22°, 9 = 24°, 9 = 40°. These an- 
gles are consistent with the values computed via Laplace's 
law directly from numerical data, that is 9 = 0.2, 9 = 
0.37, 9 = 0.69. 

Velocity profiles taken at time t — 50000 at differ- 
ent positions are shown in fig. 21 for the standard case 
H = 30, n w — 1.0. Some comments are in order. The 
velocity profile is parabolic everywhere except very near 
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Fig. 4. Velocity profile U x (y) for different cuts taken at time 
t = 40000 with the front located at at x « 210. One cut is 
taken far behind the front, x — 50, another is far ahead at 
x = 350. For these cases, approximately the same Poiseuille 
parabolic flow is found. The other two curves correspond to 
the velocities just ahead and behind the interface. In these 
case, the velocity profile is necessarily distorted in order to let 
the interface advance with an uniform velocity along y. The 
interface acts as an obstacle and the velocity shows a corre- 
sponding decrease (but not a recirculation) in the middle of 
the channel, giving rise to a two-humped profile. 
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Fig. 5. Velocity streamlines. The value of velocities are mag- 
nified by a factor 1000. The interface is located at i ~ 710. 
Near the interface the profile is distorted and a secondary flow 
appears. 



the interface. This is consistent with the assumption of a 
parabolic (Poiseuille) velocity profile. A small difference is 
present between the parabolic profile ahead and past the 
interface. This is tentatively interpreted as due to the dif- 
ferent boundary conditions applied to the fluids {n w = 1 
for the hydrophilic invading fluid 1, and n w = for fluid 
2 ahead of the front). This difference were found to dis- 
appear by setting nearer values of n w for both fluids. In 
other terms, boundary conditions are such that the fluid 
after the interface is less slipping, with a velocity at the 
wall almost recovering no-slip condition. 

In fig. [5J velocity patterns are presented. Consistently 
with fig. 4, this figure shows that the flow is one-directional 
far from the interface, confirming the assumption of a 



Poiseuille flow. Moreover, although the flow appears to 
be distorted near the interface to allow slippage, no recir- 
culation is observed at variance with other methods LBEs 
[r51fT2lfTT] . spurious currents are negligible. The spikes in 
fig. 4 reflect the existence of a hydrodynamic singularity 
near the wall. A detailed understanding of the LB dynam- 
ics in the near vicinity of this singularity remains an open 
issue for future research. 



3 Conclusions 

The present study shows that Lattice Boltzmann models 
with pseudo-potential energy interactions are capable of 
reproducing the basic features of capillary filling for binary 
fluids, as described within the Washburn approximation. 
Moreover, it has been shown that the method is able to 
reproduce the expected front dynamics for different de- 
gree of surface wettability, as well as the correct Poiseuille 
velocity profile, in the whole domain, except for a thin re- 
gion near the interface. Quantitative agreement has been 
obtained with a sufficiently thin interface, w / H < 0.3 and 
with two fluids at the same density. It would be desirable 
to extend the LB scheme in such a way to achieve larger 
density contrasts and interface widths of the order of the 
lattice spacing Ax (current values are about 5 Ax). Work 
along these lines is underway. 
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